source("setup.R")
ww1 = readRDS(fs::path("../",controls$savepoint,"ww1.rds"))
shapes = readRDS(fs::path("../",controls$savepoint,"shapes.rds"))

We continue with model A5, now extending to multiple ARAs together. We start with NUTS-2 regions.

Mittelland

We start with the Mittelland region, that includes 24 ARAs.

# select NUTS2 region
select_NUTS2 = "Mittelland"
ww_reg = ww1 %>%
  # log
  dplyr::mutate(logvl=log(vl)) %>%
  # select one ARA per NUTS-2
  dplyr::filter(NUTS2_name==select_NUTS2) %>%
  # replace zero values by 1
  dplyr::mutate(vl=if_else(vl==0, 1, vl)) %>%
  dplyr::filter(!is.na(vl)) %>%
  # create indexes for INLA
  dplyr::mutate(day1=day,
                ara1=as.numeric(as.factor(ara_n)),
                ara2=ara1)

# correspondence table
corr_reg = ww_reg %>% 
  group_by(ara_n,ara_id,ara_name,kt,pop,lab,lab_n,ara1,ara2) %>% 
  count() %>% 
  ungroup()

mw_100_desc_table(ww_reg) %>% 
  dplyr::mutate(across(everything(),as.character)) %>% 
  tidyr::gather() %>% 
  flextable::flextable(cwidth=c(4,4)) 

key

value

Number of ARAs

24

Number of laboratories

3

Number of laboratory methods

4

Number of measurements

4554

Measurements below LOQ

51

Measurements below LOD

77

First

2022-02-07

Last

2023-05-10

Median viral concentration [gc/L]

6e+04 (range: 0 to 4e+06)

Median flow [m3/day]

2e+04 (range: 4e+02 to 8e+05)

Median viral load [gc/day/100,000]

3e+12 (range: 1 to 1e+14)

mw_100_desc_table(ww_reg,ara_name) %>% 
  dplyr::select(1:7)  %>% 
  flextable::flextable(cwidth=rep(4,7)) 

ara_name

Number of ARAs

Number of laboratories

Number of laboratory methods

Number of measurements

Measurements below LOQ

Measurements below LOD

Aarwangen (Zala)

1

1

1

137

0

2

Biel

1

1

1

135

0

0

Burgdorf

1

1

1

54

0

0

Colombier (La Saunerie)

1

1

1

137

4

0

Delemont (Sede)

1

1

1

279

0

0

Fribourg (Aele)

1

1

1

329

0

8

Grenchen (Buerenamt)

1

1

1

119

0

0

Grindelwald

1

1

1

200

0

3

Interlaken

1

1

1

134

0

0

La-Chaux-de-Fonds

1

1

1

136

5

0

Laupen

1

1

1

447

35

0

Lauterbrunnen

1

1

1

143

0

1

Lyss

1

1

1

133

0

2

Mittl. Emmental

1

1

1

134

0

0

Moossee-Urtenenbach

1

1

1

136

0

20

Neuchatel

1

1

1

321

7

0

Porrentruy (Sepe)

1

1

1

228

0

2

Region Bern

1

1

1

280

0

11

Saanen

1

1

1

151

0

25

Thunersee

1

1

1

336

0

0

Vuippens (Ais)

1

1

1

134

0

0

Winznau (Zv Olten)

1

1

1

134

0

0

Worblental

1

1

1

126

0

1

Zuchwil (Soloth.-Emme)

1

1

1

191

0

2

mw_110_map_missing(ww_reg,shapes) + labs(title="Number of measurements")

ggplot(ww_reg) + geom_line(aes(x=date,y=logvl,colour=ara_name),alpha=.6) + scale_colour_discrete(guide="none") + labs(title="Viral load measurements")

Model A5.1: fixed effect

We add a geographical structure to model A5, starting by simply adding a fixed effect by ARA, assuming that the time trends of viral load across ARAs are identical. In Mittelland there is no change of methods, and there are 3 laboratories including 2 that only process one ARA each, so we leave these aspects out for now.

if(controls$rerun_models) {
  ma5.1 = INLA::inla(vl ~ 1 +
                       f(below_loq,model="iid") +
                       f(below_lod,model="iid") +
                       f(day,model="rw1", scale.model=TRUE, constr=TRUE,
                         hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
                       f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
                       f(hol,model="linear",mean.linear=0,prec.linear=.2) +
                       ara_name,
                     data = ww_reg,
                     family = "gamma",
                     control.compute = list(waic=TRUE,config=TRUE),
                     control.predictor = list(compute=TRUE,link=1))
  saveRDS(ma5.1,file=paste0("../",controls$savepoint,"ma5.1.rds"))
} else {
  ma5.1 = readRDS(file=paste0("../",controls$savepoint,"ma5.1.rds"))
}
summary(ma5.1)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", " 
##    blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
## Time used:
##     Pre = 0.71, Running = 4.95, Post = 0.186, Total = 5.84 
## Fixed effects:
##                                   mean       sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)                     13.796 2319.820  -4654.827   13.797   4682.410 13.799   0
## ara_nameBiel                     0.517    0.083      0.354    0.517      0.679  0.517   0
## ara_nameBurgdorf                 0.023    0.111     -0.195    0.023      0.241  0.023   0
## ara_nameColombier (La Saunerie)  0.368    0.083      0.206    0.368      0.531  0.368   0
## ara_nameDelemont (Sede)          0.082    0.073     -0.061    0.082      0.225  0.082   0
## ara_nameFribourg (Aele)          0.120    0.070     -0.018    0.120      0.257  0.120   0
## ara_nameGrenchen (Buerenamt)     0.300    0.086      0.131    0.300      0.470  0.300   0
## ara_nameGrindelwald              0.821    0.078      0.669    0.821      0.973  0.821   0
## ara_nameInterlaken               0.787    0.084      0.623    0.787      0.952  0.787   0
## ara_nameLa-Chaux-de-Fonds       -0.305    0.083     -0.468   -0.305     -0.142 -0.305   0
## ara_nameLaupen                   0.698    0.069      0.564    0.698      0.833  0.698   0
## ara_nameLauterbrunnen            1.448    0.083      1.286    1.448      1.610  1.448   0
## ara_nameLyss                     0.381    0.084      0.216    0.381      0.546  0.381   0
## ara_nameMittl. Emmental         -0.259    0.084     -0.423   -0.259     -0.095 -0.259   0
## ara_nameMoossee-Urtenenbach      0.125    0.084     -0.040    0.125      0.290  0.125   0
## ara_nameNeuchatel                0.677    0.071      0.538    0.677      0.816  0.677   0
## ara_namePorrentruy (Sepe)        0.062    0.075     -0.085    0.062      0.209  0.062   0
## ara_nameRegion Bern             -0.423    0.072     -0.565   -0.423     -0.282 -0.423   0
## ara_nameSaanen                   0.257    0.082      0.097    0.257      0.417  0.257   0
## ara_nameThunersee                0.186    0.070      0.049    0.186      0.323  0.186   0
## ara_nameVuippens (Ais)           0.032    0.084     -0.132    0.032      0.196  0.032   0
## ara_nameWinznau (Zv Olten)       0.138    0.084     -0.026    0.138      0.302  0.138   0
## ara_nameWorblental               0.024    0.084     -0.141    0.024      0.189  0.024   0
## ara_nameZuchwil (Soloth.-Emme)   0.337    0.078      0.184    0.337      0.490  0.337   0
## weekend                          0.007    0.030     -0.052    0.007      0.067  0.007   0
## hol                              0.075    0.078     -0.078    0.075      0.228  0.075   0
## 
## Random effects:
##   Name     Model
##     below_loq IID model
##    below_lod IID model
##    day RW1 model
## 
## Model hyperparameters:
##                                                mean    sd 0.025quant 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 2.18 0.044       2.09    2.178      2.266 2.18
## Precision for below_loq                        9.14 5.036       2.42    8.126     21.708 6.06
## Precision for below_lod                        0.00 0.000       0.00    0.000      0.000 0.00
## Precision for day                              0.33 0.045       0.25    0.326      0.428 0.32
## 
## Watanabe-Akaike information criterion (WAIC) ...: 267148.50
## Effective number of parameters .................: 542.74
## 
## Marginal log-Likelihood:  -134582.16 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
summary_exp_vl(ma5.1,pars="lab|method|hol|weekend")
##         exp(beta) 0.025quant 0.975quant
## weekend      1.01       0.95       1.07
## hol          1.08       0.93       1.26
ppp_vl_ara(ww_reg,ma5.1) 

print("Relative viral load by ARA compared to reference 'Aarwangen (Zala)':")
## [1] "Relative viral load by ARA compared to reference 'Aarwangen (Zala)':"
summary_exp_vl(ma5.1,pars="ara_n",order=TRUE)
##                                 exp(beta) 0.025quant 0.975quant
## ara_nameLauterbrunnen                4.25       3.62       5.00
## ara_nameGrindelwald                  2.27       1.95       2.65
## ara_nameInterlaken                   2.20       1.86       2.59
## ara_nameLaupen                       2.01       1.76       2.30
## ara_nameNeuchatel                    1.97       1.71       2.26
## ara_nameBiel                         1.68       1.43       1.97
## ara_nameLyss                         1.46       1.24       1.73
## ara_nameColombier (La Saunerie)      1.45       1.23       1.70
## ara_nameZuchwil (Soloth.-Emme)       1.40       1.20       1.63
## ara_nameGrenchen (Buerenamt)         1.35       1.14       1.60
## ara_nameSaanen                       1.29       1.10       1.52
## ara_nameThunersee                    1.20       1.05       1.38
## ara_nameWinznau (Zv Olten)           1.15       0.97       1.35
## ara_nameFribourg (Aele)              1.13       0.98       1.29
## ara_nameMoossee-Urtenenbach          1.13       0.96       1.34
## ara_nameDelemont (Sede)              1.09       0.94       1.25
## ara_namePorrentruy (Sepe)            1.06       0.92       1.23
## ara_nameVuippens (Ais)               1.03       0.88       1.22
## ara_nameBurgdorf                     1.02       0.82       1.27
## ara_nameWorblental                   1.02       0.87       1.21
## ara_nameMittl. Emmental              0.77       0.65       0.91
## ara_nameLa-Chaux-de-Fonds            0.74       0.63       0.87
## ara_nameRegion Bern                  0.65       0.57       0.75

We observe a large heterogeneity across ARAs, with highest relative viral load compared to reference in Lauterbrunnen, Interlaken and Grindelwald and lower relative viral load in Bern, La Chaux-de-Fonds and Emmental (note that viral load already accounts for population size).

Model A5.2: random intercept

We consider using a random intercept by ARA.

if(controls$rerun_models) {
  ma5.2 = INLA::inla(vl ~ 1 +
                       f(below_loq,model="iid") +
                       f(below_lod,model="iid") +
                       f(day,model="rw2", scale.model=TRUE, constr=TRUE,
                         hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
                       f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
                       f(hol,model="linear",mean.linear=0,prec.linear=.2) +
                       f(ara1,model="iid"),
                     data = ww_reg,
                     family = "gamma",
                     control.compute = list(waic=TRUE,config=TRUE),
                     control.predictor = list(compute=TRUE,link=1))
  saveRDS(ma5.2,file=paste0("../",controls$savepoint,"ma5.2.rds"))
} else {
  ma5.2 = readRDS(file=paste0("../",controls$savepoint,"ma5.2.rds"))
}
summary(ma5.2)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", " 
##    blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
## Time used:
##     Pre = 1.98, Running = 69.5, Post = 0.145, Total = 71.6 
## Fixed effects:
##               mean         sd  0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) 21.259 328736.466 -649479.349   21.005 649523.338 20.517   0
## weekend     -0.002      0.023      -0.047   -0.002      0.042 -0.002   0
## hol          0.103      0.066      -0.026    0.103      0.231  0.103   0
## 
## Random effects:
##   Name     Model
##     below_loq IID model
##    below_lod IID model
##    day RW2 model
##    ara1 IID model
## 
## Model hyperparameters:
##                                                mean    sd 0.025quant 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 2.04 0.011       2.02     2.04       2.06 2.04
## Precision for below_loq                        2.36 0.017       2.32     2.36       2.39 2.36
## Precision for below_lod                        0.00 0.000       0.00     0.00       0.00 0.00
## Precision for day                              0.01 0.000       0.01     0.01       0.01 0.01
## Precision for ara1                             2.48 0.020       2.44     2.48       2.52 2.47
## 
## Watanabe-Akaike information criterion (WAIC) ...: 267513.60
## Effective number of parameters .................: 557.74
## 
## Marginal log-Likelihood:  -136432.44 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
summary_exp_vl(ma5.2,pars="lab|method|hol|weekend")
##         exp(beta) 0.025quant 0.975quant
## weekend      1.00       0.95       1.04
## hol          1.11       0.97       1.26
ppp_vl_ara(ww_reg,ma5.2) 

print("Relative viral load by ARA compared to average:")
## [1] "Relative viral load by ARA compared to average:"
ma5.2$summary.random$ara1 %>% 
  as_tibble() %>% 
  dplyr::transmute(ara1=ID,
                   `exp(beta)`=round(exp(mean),2),
                   `0.025quant`=round(exp(`0.025quant`),2),
                   `0.975quant`=format(round(exp(`0.975quant`),2),scientific=FALSE)) %>% 
  dplyr::left_join(select(corr_reg,ara1,ara_name),by = join_by(ara1)) %>% 
  dplyr::arrange(-`exp(beta)`) %>% 
  select(-ara1) %>% 
  column_to_rownames(var="ara_name")
##                         exp(beta) 0.025quant 0.975quant
## Lauterbrunnen                3.20       2.42       4.23
## Grindelwald                  1.75       1.33       2.30
## Interlaken                   1.67       1.26       2.20
## Laupen                       1.50       1.15       1.95
## Neuchatel                    1.49       1.14       1.94
## Biel                         1.26       0.96       1.67
## Lyss                         1.14       0.86       1.50
## Colombier (La Saunerie)      1.10       0.83       1.45
## Zuchwil (Soloth.-Emme)       1.08       0.82       1.41
## Grenchen (Buerenamt)         1.02       0.77       1.36
## Saanen                       0.98       0.74       1.30
## Thunersee                    0.91       0.70       1.19
## Winznau (Zv Olten)           0.89       0.67       1.18
## Delemont (Sede)              0.88       0.67       1.15
## Fribourg (Aele)              0.87       0.67       1.14
## Moossee-Urtenenbach          0.86       0.65       1.14
## Porrentruy (Sepe)            0.82       0.62       1.07
## Vuippens (Ais)               0.79       0.60       1.05
## Worblental                   0.79       0.60       1.05
## Burgdorf                     0.78       0.57       1.06
## Aarwangen (Zala)             0.76       0.57       1.00
## Mittl. Emmental              0.60       0.45       0.79
## La-Chaux-de-Fonds            0.58       0.44       0.77
## Region Bern                  0.51       0.39       0.67

Results are similar in terms of rankings. Of course the reference is now the between-ARA mean instead of one arbitrary ARA, which is more easily interpreted. As expected ARA-level intercepts are pulled towards the mean.

Model A5.3: space-time interaction

We now allow different time trends across ARAs.

if(controls$rerun_models) {
  ma5.3 = INLA::inla(vl ~ 1 +
                       f(below_loq,model="iid") +
                       f(below_lod,model="iid") +
                       f(day,model="rw2", scale.model=TRUE, constr=TRUE,
                         hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
                       f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
                       f(hol,model="linear",mean.linear=0,prec.linear=.2) +
                       f(ara1,model="iid") +
                       f(day1,model="rw1", scale.model=TRUE, constr=TRUE,
                         group=ara2, control.group=list(model="iid"),
                         hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))),
                     data = ww_reg,
                     family = "gamma",
                     control.compute = list(waic=TRUE,config=TRUE),
                     control.predictor = list(compute=TRUE,link=1))
  saveRDS(ma5.3,file=paste0("../",controls$savepoint,"ma5.3.rds"))
} else {
  ma5.3 = readRDS(file=paste0("../",controls$savepoint,"ma5.3.rds"))
}
summary(ma5.3)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", " 
##    blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
## Time used:
##     Pre = 0.798, Running = 110, Post = 0.588, Total = 111 
## Fixed effects:
##               mean     sd 0.025quant 0.5quant 0.975quant   mode   kld
## (Intercept) 13.899 10.975     -9.345   13.899     37.137 13.898 0.005
## weekend     -0.001  0.020     -0.040   -0.001      0.037 -0.001 0.000
## hol          0.104  0.058     -0.009    0.104      0.217  0.104 0.000
## 
## Random effects:
##   Name     Model
##     below_loq IID model
##    below_lod IID model
##    day RW2 model
##    ara1 IID model
##    day1 RW1 model
## 
## Model hyperparameters:
##                                                  mean    sd 0.025quant 0.5quant 0.975quant   mode
## Precision parameter for the Gamma observations  2.908 0.072      2.770    2.907      3.054  2.903
## Precision for below_loq                         0.968 1.195      0.038    0.528      4.969  0.119
## Precision for below_lod                         0.005 0.002      0.002    0.005      0.011  0.004
## Precision for day                               0.010 0.001      0.009    0.010      0.013  0.010
## Precision for ara1                             11.751 2.064      7.697   11.585     15.971 11.784
## Precision for day1                              1.166 0.119      0.946    1.161      1.414  1.153
## 
## Watanabe-Akaike information criterion (WAIC) ...: 265619.78
## Effective number of parameters .................: 634.31
## 
## Marginal log-Likelihood:  -159294.52 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
summary_exp_vl(ma5.3,pars="lab|method|hol|weekend")
##         exp(beta) 0.025quant 0.975quant
## weekend      1.00       0.96       1.04
## hol          1.11       0.99       1.24
ppp_vl_ara(ww_reg,ma5.3)

print("Relative viral load by ARA compared to average:")
## [1] "Relative viral load by ARA compared to average:"
ma5.3$summary.random$ara1 %>% 
  as_tibble() %>% 
  dplyr::transmute(ara1=ID,
                   `exp(beta)`=round(exp(mean),2),
                   `0.025quant`=round(exp(`0.025quant`),2),
                   `0.975quant`=format(round(exp(`0.975quant`),2),scientific=FALSE)) %>% 
  dplyr::left_join(select(corr_reg,ara1,ara_name),by = join_by(ara1)) %>% 
  dplyr::arrange(-`exp(beta)`) %>% 
  select(-ara1) %>% 
  column_to_rownames(var="ara_name")
##                         exp(beta) 0.025quant 0.975quant
## Laupen                       1.63       1.38       1.92
## Lauterbrunnen                1.52       1.04       2.30
## Grindelwald                  1.48       1.25       1.76
## Neuchatel                    1.43       1.21       1.70
## Biel                         1.24       0.86       1.83
## Interlaken                   1.21       0.82       1.80
## Zuchwil (Soloth.-Emme)       1.14       0.96       1.35
## Lyss                         1.10       0.75       1.62
## Saanen                       1.05       0.73       1.51
## Colombier (La Saunerie)      1.04       0.71       1.52
## Grenchen (Buerenamt)         1.00       0.68       1.46
## Winznau (Zv Olten)           0.96       0.66       1.40
## Thunersee                    0.95       0.81       1.13
## Vuippens (Ais)               0.93       0.64       1.34
## Burgdorf                     0.92       0.56       1.51
## Moossee-Urtenenbach          0.91       0.62       1.34
## Worblental                   0.90       0.62       1.31
## Fribourg (Aele)              0.88       0.75       1.05
## Aarwangen (Zala)             0.87       0.59       1.28
## Porrentruy (Sepe)            0.85       0.72       1.01
## Delemont (Sede)              0.82       0.70       0.97
## La-Chaux-de-Fonds            0.72       0.48       1.05
## Mittl. Emmental              0.70       0.47       1.02
## Region Bern                  0.52       0.44       0.62
ndays = length(unique(ww_reg$day))
nara = length(unique(ww_reg$ara1))
ma5.3$summary.random$day1 %>% 
  bind_cols(day=rep(0:(ndays-1),nara),
            ara1=rep(1:nara,each=ndays)) %>% 
  left_join(corr_reg,by = join_by(ara1)) %>% 
  ggplot() +
  geom_hline(yintercept=1,linetype=2,alpha=.5) +
  geom_ribbon(aes(x=day,ymin=exp(`0.025quant`),ymax=exp(`0.975quant`),fill=ara_name),alpha=.5) +
  geom_line(aes(x=day,y=exp(mean),colour=ara_name)) +
  facet_wrap(~ara_name) +
  scale_colour_discrete(guide="none") +
  scale_fill_discrete(guide="none") +
  scale_y_continuous(trans="log",breaks = c(.1,1,10)) +
  coord_cartesian(ylim=c(.05,20)) +
  labs(title="Deviations from average time trend by ARA",x="Day",y="Relative viral load by ARA") 

This brings a clear improvement in WAIC. We observe similar average between-ARA heterogeneity, with highest relative viral loads in Laupen, Lauterbrunnen, Grindelwald and Neuchâtel (four ARAs where the credible intervals are the highest) and lowest relative viral loads in Delemont, La-Chaux-de-Fonds, Emmental and Bern (four ARAs where the credible intervals are the lowest). While time trends are generally aligned, we observe deviations from the regional average time trend in some ARAs such as Colombier, Delemont and La-Chaux-de-Fonds (comparatively lower viral loads during Summer 2022), and in Lauterbrunnen, Grindelwald and Interlaken (comparatively higher during Summer 2022).

Model A5.4: spatial structure

We now consider the neighboring structure between ARAs.

# setup neighboring matrix
shapes_reg = shapes$ara_shp %>% 
  dplyr::filter(ara_id %in% ww_reg$ara_id) %>% 
  left_join(corr_reg,by = join_by(ara_id)) %>% 
  dplyr::arrange(ara1)
sf_use_s2(FALSE)
graph_reg = spdep::poly2nb(shapes_reg)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
path_graph = paste0("../",controls$savepoint,"W_reg_",select_NUTS2,".adj")
nb2INLA(path_graph, graph_reg)
if(controls$rerun_models) {
  ma5.4 = INLA::inla(vl ~ 1 +
                       f(below_loq,model="iid") +
                       f(below_lod,model="iid") +
                       f(day,model="rw2", scale.model=TRUE, constr=TRUE,
                         hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
                       f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
                       f(hol,model="linear",mean.linear=0,prec.linear=.2) +
                       f(ara1,model="bym2",
                         graph=path_graph,
                         scale.model = TRUE, constr = TRUE, 
                         hyper = list(theta1 = list("PCprior", c(1, 0.01)),  # Pr(sd<1) = 0.01, unlikely to have rr>3just based on the spatial confounding
                                      theta2 = list("PCprior", c(0.5, 0.5)))  # Pr(phi<0.5)=0.5, we state that we believe that the unmeasured spatial confounding is driven 50% from the structured and 50% from the unstructured random effect
                       ) +
                       f(day1,model="rw1", scale.model=TRUE, constr=TRUE,
                         group=ara2, control.group=list(model="iid"),
                         hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))),
                     data = ww_reg,
                     family = "gamma",
                     control.compute = list(waic=TRUE,config=TRUE),
                     control.predictor = list(compute=TRUE,link=1))
  saveRDS(ma5.4,file=paste0("../",controls$savepoint,"ma5.4.rds"))
} else {
  ma5.4 = readRDS(file=paste0("../",controls$savepoint,"ma5.4.rds"))
}
summary(ma5.4)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", " 
##    blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
## Time used:
##     Pre = 8.61, Running = 142, Post = 1.54, Total = 152 
## Fixed effects:
##               mean     sd 0.025quant 0.5quant 0.975quant   mode   kld
## (Intercept) 13.908 12.669    -11.510   13.908     39.319 13.908 0.002
## weekend     -0.002  0.020     -0.040   -0.002      0.037 -0.002 0.000
## hol          0.104  0.058     -0.009    0.104      0.217  0.104 0.000
## 
## Random effects:
##   Name     Model
##     below_loq IID model
##    below_lod IID model
##    day RW2 model
##    ara1 BYM2 model
##    day1 RW1 model
## 
## Model hyperparameters:
##                                                  mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision parameter for the Gamma observations  2.913 0.073      2.767    2.914      3.056 2.919
## Precision for below_loq                         3.816 3.145      0.474    2.971     12.057 1.372
## Precision for below_lod                         0.005 0.008      0.000    0.003      0.025 0.001
## Precision for day                               0.011 0.002      0.007    0.010      0.016 0.010
## Precision for ara1                             10.756 4.008      5.168   10.005     20.694 8.673
## Phi for ara1                                    0.183 0.182      0.008    0.118      0.691 0.020
## Precision for day1                              1.202 0.148      0.949    1.188      1.533 1.154
## 
## Watanabe-Akaike information criterion (WAIC) ...: 265628.62
## Effective number of parameters .................: 638.86
## 
## Marginal log-Likelihood:  -159271.80 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
summary_exp_vl(ma5.4,pars="lab|method|hol|weekend")
##         exp(beta) 0.025quant 0.975quant
## weekend      1.00       0.96       1.04
## hol          1.11       0.99       1.24
ppp_vl_ara(ww_reg,ma5.4) 

print("Relative viral load by ARA compared to average:")
## [1] "Relative viral load by ARA compared to average:"
ma5.4$summary.random$ara1 %>% 
  as_tibble() %>% 
  dplyr::transmute(ara1=ID,
                   `exp(beta)`=round(exp(mean),2),
                   `0.025quant`=round(exp(`0.025quant`),2),
                   `0.975quant`=format(round(exp(`0.975quant`),2),scientific=FALSE)) %>% 
  dplyr::left_join(select(corr_reg,ara1,ara_name),by = join_by(ara1)) %>% 
  dplyr::arrange(-`exp(beta)`) %>% 
  dplyr::select(-ara1) %>% 
  dplyr::filter(!is.na(ara_name)) %>% 
  column_to_rownames(var="ara_name")
##                         exp(beta) 0.025quant 0.975quant
## Lauterbrunnen                1.64       1.10       2.47
## Laupen                       1.62       1.38       1.92
## Grindelwald                  1.48       1.25       1.76
## Neuchatel                    1.43       1.22       1.69
## Interlaken                   1.30       0.87       1.95
## Biel                         1.28       0.85       1.94
## Zuchwil (Soloth.-Emme)       1.13       0.96       1.34
## Lyss                         1.07       0.73       1.58
## Saanen                       1.05       0.72       1.53
## Grenchen (Buerenamt)         1.00       0.68       1.48
## Colombier (La Saunerie)      0.99       0.67       1.46
## Thunersee                    0.95       0.81       1.12
## Winznau (Zv Olten)           0.95       0.64       1.41
## Vuippens (Ais)               0.91       0.62       1.34
## Fribourg (Aele)              0.88       0.75       1.04
## Moossee-Urtenenbach          0.86       0.59       1.27
## Burgdorf                     0.85       0.51       1.41
## Porrentruy (Sepe)            0.85       0.72       1.00
## Worblental                   0.84       0.58       1.24
## Aarwangen (Zala)             0.83       0.56       1.24
## Delemont (Sede)              0.82       0.69       0.97
## La-Chaux-de-Fonds            0.70       0.46       1.03
## Mittl. Emmental              0.67       0.45       0.98
## Region Bern                  0.52       0.44       0.62
ndays = length(unique(ww_reg$day))
nara = length(unique(ww_reg$ara1))
ma5.4$summary.random$day1 %>% 
  bind_cols(day=rep(0:(ndays-1),nara),
            ara1=rep(1:nara,each=ndays)) %>% 
  left_join(corr_reg,by = join_by(ara1)) %>% 
  ggplot() +
  geom_hline(yintercept=1,linetype=2,alpha=.5) +
  geom_ribbon(aes(x=day,ymin=exp(`0.025quant`),ymax=exp(`0.975quant`),fill=ara_name),alpha=.5) +
  geom_line(aes(x=day,y=exp(mean),colour=ara_name)) +
  facet_wrap(~ara_name) +
  scale_colour_discrete(guide="none") +
  scale_fill_discrete(guide="none") +
  scale_y_continuous(trans="log",breaks = c(.1,1,10)) +
  coord_cartesian(ylim=c(.05,20)) +
  labs(title="Deviations from average time trend by ARA",x="Day",y="Relative viral load by ARA") 

Results are very similar to model A5.3. We find that the spatial structure explains between 0 and 70% of the spatial heterogeneity.